install.packages("DESeq2")
BiocManager::install("DESeq2")
library("DESeq2")
BiocManager::install("limma")
library("limma")
BiocManager::install("edgeR")
library("edgeR")
library(Biobase)
library(DESeq2)
# 加载所需包
library(DESeq2)



count<-read.csv("count.csv", header = TRUE, stringsAsFactors = FALSE, sep = " ", row.names = 1)
clin<-read.csv("clin_inf.csv", header = TRUE, stringsAsFactors = FALSE, sep = " ", row.names = 1)
ID<-read.csv("ID_annoation.csv", header = TRUE, stringsAsFactors = FALSE, sep = " ", row.names = 1)
exp<-read.csv("exp_inf.csv", header = TRUE, stringsAsFactors = FALSE, sep = " ", row.names = 1)

meta<-merge(exp,clin,by.x=0,by.y=0,all=FALSE)
rownames(meta)<-meta$Row.names
meta<-meta[,-1]

library(Biobase)
x<-colnames(count)
x<-gsub(".$","",x)
x<-gsub("X(...$)","PTB\\1",x)
x<-gsub("X(.*)","XYA\\1",x)
colnames(count)<-x

y<-rownames(ID)
y<-gsub("\\..*","",y)
a<-!duplicated(y)
ID<-ID[a,]
rownames(ID)<-y[a]
ID<-ID[match(rownames(count),rownames(ID)),]

count<-count[,match(rownames(meta),colnames(count))]
#创建Exprssionset对象
eset<-ExpressionSet(assayData = as.matrix(count),
                    phenoData = new("AnnotatedDataFrame",data=meta),
                    featureData=new("AnnotatedDataFrame",data=ID))
eset2 <- eset[1:10,1:10]
# 筛选出 Treatment 为 "Treatment" 的样本
treated_samples <- eset[, eset$批次 == "batch1"]
# 筛选出位于 "chr1" 的基因
chrX_genes <- eset[fData(eset)$seqnames == "chr1", ]
meta$age_group <- ifelse(meta$年龄 < 65, "middle_age", "elderly")
coldata <- data.frame(
  condition = meta$age_group,
  row.names = rownames(meta)
)
count<-round(count)
dds <- DESeqDataSetFromMatrix(
  countData = count,
  colData = coldata,
  design = ~ condition
)
dds <- DESeq(dds)
res<-results(dds)
# 处理 DESeq2 结果
# 处理 DESeq2 结果
res_df <- as.data.frame(res)
# 过滤掉 padj 或 log2FoldChange 为 NA 的行
res_df <- res_df[!is.na(res_df$padj) &!is.na(res_df$log2FoldChange), ]
res_df$sig <- "NS"
res_df$sig[res_df$padj < 0.05 & abs(res_df$log2FoldChange) > 1] <- ifelse(res_df$log2FoldChange[res_df$padj < 0.05 & abs(res_df$log2FoldChange) > 1] > 0, "Up", "Down")
df1<-read.csv(file = "class2/clin_inf.csv",sep = " ",header=TRUE,row.names = 1)
df2<-read.csv(file = "class2/count.csv",sep = " ",header=TRUE,row.names = 1)
df3<-read.csv(file = "class2/exp_inf.csv",sep = " ",header=TRUE,row.names = 1)
df4<-read.csv(file = "class2/ID_annoation.csv",sep = " ",header=TRUE)
col1=colnames(df2)
col2<- gsub(".$","\\1",col1)
col3<-gsub("X(...$)","PTB\\1",col2)
col4<-gsub("X","XYA",col3)
colnames(df2)<-col4
row2<-df4[,1]
row3<-gsub("\\..*", "", row2) 
df4[,1]<-row3
df4<-df4[!duplicated(df4$gene_id),]
rownames(df4) <- df4[, 1]
df4 <- df4[, -1]
df5<-merge(df1,df3,by.x=0,by.y=0)
rownames(df5) <- df5[, 1]
df5<- df5[, -1]
df5 <- df5[match(colnames(df2),rownames(df5)),]
df4<-df4[match(rownames(df2),rownames(df4)),]
eset <- ExpressionSet(assayData = as.matrix(df2),
                      phenoData = new("AnnotatedDataFrame", data = df5),
                      featureData = new("AnnotatedDataFrame", data = df4))

df5$age_group <- ifelse(df5$年龄 < 65, "middle_age", "elderly")
coldata <- data.frame(
  condition = df5$age_group,
  row.names = rownames(df5)
)
df2<-round(df2)
dds <- DESeqDataSetFromMatrix(
  countData = df2,
  colData = coldata,
  design = ~ condition
)
dds <- DESeq(dds)
res<-results(dds)
# 处理 DESeq2 结果
# 处理 DESeq2 结果
res_df <- as.data.frame(res)
# 过滤掉 padj 或 log2FoldChange 为 NA 的行
res_df <- res_df[!is.na(res_df$padj) &!is.na(res_df$log2FoldChange), ]
res_df$sig <- "NS"
res_df$sig[res_df$padj < 0.05 & abs(res_df$log2FoldChange) > 1] <- ifelse(res_df$log2FoldChange[res_df$padj < 0.05 & abs(res_df$log2FoldChange) > 1] > 0, "Up", "Down")

# 绘制火山图
ggplot<-ggplot(res_df, aes(x = log2FoldChange, y = -log10(padj))) +
  geom_point(aes(color = sig)) +
  scale_color_manual(values = c("NS" = "gray", "Up" = "red", "Down" = "blue")) +
  labs(x = "Log2 Fold Change", y = "-Log10 Adjusted p-value", title = "Volcano Plot") +
  theme_minimal() +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed") +
  geom_vline(xintercept = c(-1, 1), linetype = "dashed")
ggsave("火山图.png", width = 8, height = 6)
counts_matrix <- counts(dds, normalized = FALSE)
# 筛选差异表达基因的行名
diff_genes <- rownames(res_df[res_df$sig != "NS", ])
# 从表达矩阵中提取差异表达基因的表达数据
diff_expr_matrix <- counts_matrix[diff_genes, ]
install.packages("pheatmap")
library(pheatmap)
pheatmap(diff_expr_matrix, 
         scale = "row", 
         cluster_rows = TRUE, 
         cluster_cols = TRUE,
         annotation_col = coldata,
         show_rownames = FALSE, 
         show_colnames = TRUE,
         color = colorRampPalette(c("blue", "white", "red"))(50),
         filename = "热图.png")
